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Abstract 

The Stieltjes constants 7 „ appear in the coefficients in the Laurent expansion 
of the Riemann zeta function C(s) about the simple pole s = 1. We present an 
asymptotic expansion for 7 „ as n —>■ oo based on the approach described by Knessl 
and Coffey [Math. Comput. 80 (2011) 379-386]. A truncated form of this expansion 
with explicit coefficients is also given. Numerical results are presented that illustrate 
the accuracy achievable with our expansion. 
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1. Introduction 

The Stieltjes constants 7 „ appear in the coefficients in the Laurent expansion of the 
Riemann zeta function ^(s) about the point s = 1 given by 

C(<>) = -r + 

“ - 1 „„ n! 

where 70 = 0.577216... is the well-known Euler-Mascheroni constant. Some historical 
notes and numerical values of 7 ^ for n < 20 are given in [3]. Recent high-precision 
evaluations of 7 ^ based on numerical integration have been described in [5, 8 ]. In [5], 
Keiper lists various 7 ^ up to n = 150, whereas in [ 8 ], Kreminski has computed values 
to several thousand digits for n < 10 ^ and for further selected values (accurate to 10 ^ 
digits) up to n = 5 X 10^. All values up to n = 10^ have been computed by Johansson in 
[4] to about 10^ digits. 

Upper bounds for | 7 „| in the form 

i7„i<{3+(-r}LlM, 

have been obtained by Berndt [1] with = 1, and by Zhang and Williams [13] with 
A„ = ( 2 /n)"’ 7 r“m(n -|- |) ~ y/2{2leY for n — 00 . On the other hand, Matsuoka [9] has 
shown that 

l7n| < 10 -V^°si°s" (n> 10 ). 
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However, all these bounds grossly overestimate the growth of for large values of n. 
An asymptotic approximation for 7 ^ has recently been given by Knessl and Coffey [ 6 ] in 
the form 

B^nA 

7n ~— 7 ^ COS (na + 6 ) (n—)-+oo), (1.1) 

yjn 

where A, B, a and b are functions that depend weakly on n; see Section 2 for the definition 
of these quantities. Knessl and Coffey have verified numerically that for n < 3.5 x 10^ 
the above formula accounts for the asymptotic growth and oscillatory pattern of 7 ^, with 
the exception of n = 137 where the cosine factor in (1.1) becomes very small. 

The aim in this note is to extend the analysis in [ 6 ] to generate an asymptotic expan¬ 
sion for 7 „ as n — 7 - -|-oo. The coefficients in this expansion are determined numerically 
by application of Wojdylo’s formulation [14] for the coefficients in the expansion of a 
Laplace-type integral. An explicit evaluation of the coefficients is obtained in the case 
of the expansion truncated after three terms. This approximation is extended to the 
more general Stieltjes constants 7 n(cr) appearing in the Laurent expansion of the Hur- 
witz zeta function C('S,a)- Numerical results are presented in Section 3 to demonstrate 
the accuracy of our expansion compared to that in ( 1 . 1 ). 


2. Asymptotic expansion for 7 ^ 

We start with the integral representation for n > 1 given in [13] 

log"^“^ X 

7n = / Bi{x - [x]) - 2 — “ log dx, 

J \ X 

where Ri(x — [a;]) = — is the first periodic Bernoulli polynomial. With the 

change of variable t = log x, we obtain [6, Eq. (2.3)] 

7„ = —A I ^ — f exp [27rffce* -|- n log t — t] f -1^ dt\. 

TTK Jo \t / 

Following the approach used in [6], we define 

A* C ^ / 

ipkit) = V’fc(L n) = - — e* - log t, f{t) = f{t; n) = 


( 2 . 1 ) 


and write 


'y ^ ’^k: Jk ■ — 


fc=l 


n 

Tik 


^-nipkd) 


f{t) dt. 


( 2 . 2 ) 


We employ the method of steepest descents to estimate the integrals Jk for large 
values of n. Saddle points of the exponential factor occur at the zeros of ipkit) = 0; that 
is, they satisfy 

te^ = 7 ^. (2.3) 


27rk 


There is an infinite string of saddle points, which is approximately parallel to the imagi¬ 
nary f-axis, given by [ 6 ] 

tm = log - loglog n + {2m + ^)7ri -f 

2'Kk \ log n J 
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Figure 1: Paths of steepest descent and ascent through the saddles to and ti when n = 100 and 
k = 1. The steepest paths through the saddle t_i (not shown) in the lower half-plane are similar 
to those through ti. The arrows indicate the direction of integration. 


for m = 0, ±1, ±2,... and large n. For fixed k and m, the value of 3ft '4>k{tm) is then 

-3ftV’fc(tm) = log log n- (1 -h log (27rA: log n)) -h 0((log n)“^) 

log n 


as n —)• oo, where the dependence on m is contained in the order term. This shows that 
the heights of the saddles corresponding to A: > 2 are exponentially smaller as n —)• oo 
than the saddle with k = 1, so that to within exponentially small correction terms we 
may neglect the contribution in (2.2) arising from k values corresponding to A: > 2; but 
see the discussion in Section 3. From hereon, we shall drop the subscript k and write 

V’i(A) = V’(A)- 

Typical paths of steepest descent and ascent through the saddles to and ti are shown 
in Fig. 1. Steepest descent and ascent paths terminate at infinity in the right-half plane 
in the directions Q{t) = {2j + i)7r and $>(t) = {2j + |)7r (j = 0, ±1, ±2,...), respectively. 
The steepest descent paths through to and ti emanate from the origin and pass to infinity 
in the directions $^(t) = ^vr and |7r, respectively. Similarly, the steepest descent path 
through t-i (not shown) emanates from the origin and passes to infinity in the direction 
A(t) = — Ivr. The integration path in (2.2) can then be deformed to pass through the 
saddle to as shown in Fig. 1. 

Application of the method of steepest descents (see, for example, [10, p. 127] and [11, 
p. 14]) then yields 


^ ^g-n^(to)-Io / ^ “ C2s(|). 

TT to(-'0"(Ao))^/^ \ ny^n*+V2’ 


(2.4) 


where (a)* = r(a -|- s)/r(a) is the Pochhammer symbol and cq = 1- The normalised 
coefficients C 2 s can be obtained by an inversion process and are listed for s < 4 in [2, 
p. 119] and for s < 2 in [11, p. 13]; see below. Alternatively, they can be obtained by an 
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expansion process to yield Wojdylo’s formula [14] given by 


_ -s/2 l^s-k + \)j 


E h's—k 

/3q ^ 

fc=0 j=0 


j! Oq 


j > 


(2.5) 


see also [12, p. 25]. Here B^j = 02 , • • •, CKfc-j+i) are the partial ordinary Bell 

polynomials generated by the recursion^ 


k-j+l 

^kj — ^ ^ ^rBk—rJ — l‘> 

r=l 


BkO — ^/cO? 


where 6mn is the Kronecker symbol, and the coefficients and j3r appear in the expan¬ 
sions 

00 00 

\r+2 Q u 


ij{t) - Ipito) = f{t) = '^Pr(t-toY (2.6) 

r=0 r=0 

valid in a neighbourhood of the saddle t = to- 

Following [6], we put to = u + iv, where u, v are real, and write —ipito) = A + ia, 
where 

A : = §f?(log to - 1/to) = ^ log(M^ + u^) - 


u 


a : = $^(log to — 1/to) = arctanf — ) + 

\uj 

We have V’^^(^o) = (1 + ^o)/^o accordingly define^ 


+ t;2 ' 

V 




(2.7) 


B := 2V^ 


to 


b := ^TT — V — arctan 


\/l + to 

A simple calculation using (2.3) with k = 1 shows that 


1 + tt 


u 

tanu = —, 
V 


27r|to| 


n 


( 2 . 8 ) 


(2.9) 


Then, from (2.2) with k = 1, (2.4) and the second relation in (2.9), we find upon incor¬ 
porating the factor 1 — to/n into the asymptotic series that 


Be 


where 


7n 


C2s — C2s 


nA ( oo , /In . 

n 1 ^ I 


s=0 


2to 


■C2s-2 (s > 1). 


2s - 1 

If we now introduce the real and imaginary parts of the coefficients C 2 s by 
C2s + iDs (s ^ 1); Co = 1, Do = 0, 


( 2 . 10 ) 

( 2 . 11 ) 


where we recall that Cg and Dg contain an n-dependence, then we have the expansion of 
7 n given by the following theorem. 

'^For example, this generates the values B41 = 04, B42 = 03 -I- 20103, B43 = 80^02 and B44 = af. 

^In [6], the quantity I-jt — v appearing in the definition of b is written as arctan (v/u) by virtue of the 
first relation in (2.9). 
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Theorem 1. Let the quantities A, B, a and h, and the coefficients Cg, Dg, be as defined 
in (2.1), (2.8) and (2.11). Then, neglecting exponentially smaller terms, we have 


In ~ 



CXD 

— sin {na + b) ^ 

S=1 


i 


( 2 . 12 ) 


as n —)• 00 . 


We note that to leading order A ~ log log n and B ~ (Svrlog for large n. 

A simpler form of the expansion (2.12) can be given by truncating the above series 
at s = 2 and use of the form of the normalised coefficients C 2 s in (2-4) expressed in the 
form 

« = (2^,4„))2 {§A - ffsA + KJ'I'l - 'I'4)F2 - f (>1-1 - >l-3>l-4 + 

+f (nl'l - K'l'a - 5'I'4)'I'4 + il-al-s - ^4-6)} 

where, for brevity, we have defined 


:= 


V{to) 


(m > 3), 


F — 

^ m • — 


fito) 


(m > 1); 


see [2, p. 119], [11, pp. 13-14]. 

From (2.1) and (2.10), use of Mathematica shows that 


^2 — 


p2(to) (4 + 3to)tg 

12(1+ to)^ n(l + to)^ 


+ 0{n ^), 


C 4 = 


P4(to) 


864(1 + to)‘^ 


+ 0{n 




where 

P 2 {to) = 2 - 18to - 20^2 - Ml + 24, 

pfito) = 4 - 72to - 3324 - 8028t;) - 19644t^ - 20280t^ - 9911t|] - 1884t^ + 4t^. 


Then we obtain the following result. 


Theorem 2. Let the quantities A, B, a and b be as defined in (2.7) and (2.8). Then, 
with 


Cl + idi 


p2{to) 

24(1+ to4' 


C2 + id2 


pijto) (4 + 3to)4 

1152(1 +to)® 2(l + to)2 ’ 


where Cg, dg (s = 1,2) are real (and independent of n) and to is the saddle point given 
by the principal solution of (2.3) with k = 1, we have the asymptotic approximation 


Be 


nA 


In 


n 


I cos (na + h)(l + — + - sin (na + b) (— + 1 (2-13) 

L \ n J \n J) 


as n —)• 00 . 
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We remark that the expansion of the integrals for fixed k > 2 follows the same 
procedure. If we still refer to the real and imaginary parts of the contributory saddle to 
(when k > 2) as u and v, the second relation in (2.9) is now replaced by = 27rA:|to|/'R- 
It then follows that the form of the expansion for —'^Jk is given by (2.12), provided the 
quantities A, B, a and b, and the coefficients Cg, Dg, are interpreted as corresponding to 
the saddle to with the /c-value under consideration. 


3. Numerical results and concluding remarks 

We discuss numerical computations carried out using the expansions given in Theorems 
1 and 2. For a given value of n the saddle to is computed from (2.3) with k = 1 to the 
desired accuracy. Mathematica is used to determine the coefficients and jSr in (2.6) 
for 0 < r < 2so, where in the present computations sq = 6. The coefficients Cg and Dg 
can then be calculated for 0 < s < sq from (2.5), (2.10) and (2.11). 

We display the computed values of Cg and Dg for two values of n in Table 1. We repeat 
that these coefficients contain an n-dependence and so have to be computed for each value 
of n chosen. In Table 2, the values of the absolute relative error in the computation of 
from the expansion (2.12) are presented as a function of the truncation index s for 
several values of n. 


Table 1 ; Values of the coefficients Cs and Ds (to 10 dp) for 1 < s < 6 for two values of n. 


s 

n = 

Cg 

100 

Dg 

II 

e 

1000 

Dg 

1 

-0.3158578918 

-h0.1626819326 

-0.0885061806 

-h0.1958085240 

2 

-2.9096870797 

-2.1947177121 

-6.5840165991 

-2.6459812815 

3 

-0.3804847598 

-3.3953890569 

-9.4682639154 

-10.09635962642 

4 

-hi.4820479884 

-0.1130053628 

-1.3074432243 

-11.31040992292 

5 

-0.2630549338 

-hO.9253656779 

-h4.9469591967 

-1.67819725309 

6 

-0.3783700609 

-0.3119889058 

-h0.8180579543 

-h3.98701271605 


The case n = 137 has been included in Table 2 since this corresponds to the factor 
cos(na + 6) possessing the very small value ~ 1.69881 x 10“^. The leading term approxi¬ 
mation in (1.1), and (2.12) (with s = 0), yields an incorrect sign, namely -1-3.89874 x 10^^ 
when 7137 = —7.99522199 ... x 10^^. According to [4], this is the only instance for n < 10® 
when the leading approximation produces the wrong sign. It is seen that inclusion of 
the higher order correction terms with s < 6 yields a relative error of order 10“^*^ in this 
case. When n = 10®, [4] gives the value 

'yiooooo = 1.99192730631254109565822724315... x 10®®^^^ 

The expansion (2.12) for this value of n with truncation index s = 6 is found to yield 
a relative error of order 10“®^; that is, the expansion correctly reproduces all the digits 
displayed above. 
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Table 2: Values of the absolute relative error in the computation of 7 n from (2.12) as a function of the 
truncation index s for different n. 


s 

n = 75 

n = 100 

n = 137 

n = 1000 

0 

1.759 X 10-3 

1.412 X 10-3 

— 

1.597 X 10-^ 

1 

6.503 X 10-^ 

3.226 X 10-^ 

2.701 X 10-1 

2.649 X 10-3 

2 

1.244 X 10-3 

4.472 X 10-3 

8.775 X 10-2 

4.125 X 10-9 

3 

3.063 X 10-^ 

9.370 X 10-3 

3.811 X 10-3 

7.711 X 10-11 

4 

2.535 X 10-9 

7.850 X 10-1° 

2.183 X 10-3 

2.026 X 10-13 

5 

5.101 X 10-1° 

9.022 X 10-11 

1.248 X 10-3 

6.157 X 10-13 

6 

1.850 X 10-13 

1.982 X 10-12 

9.415 X 10-1° 

2.743 X 10-13 


Table 3: Values of the absolute relative error in the computation of 7 „ from (2.12) with fc = 1 and 
fc < 2 as a function of the truncation index s for n — 25. 


s 

k = 1 

k<2 

0 

1.051 X 10-2 

1.052 X 10-2 

1 

2.909 X 10-3 

2.894 X 10-3 

2 

2.608 X 10-^ 

2.460 X 10-^ 

3 

2.390 X 10-3 

1.723 X 10-3 

4 

1.518 X 10-3 

3.412 X 10-'^ 

5 

1.495 X 10-3 

1.160 X 10-'^ 

6 

1.482 X 10-3 

1.189 X 10-3 


For the smallest value n = 75 presented in Table 2, it is found numerically that the 
contribution to (2.2) corresponding to /c = 2 is about 11 orders of magnitude smaller than 
the dominant term with k = 1. For the larger n values, this contribution is even smaller 
and the terms with k >2 can be safely neglected. However, for smaller n this is no longer 
the case and a meaningful approximation has to take into account the contribution from 
other k > 2 values. 

In Table 3, we illustrate this situation by taking n = 25. The second column shows 
the absolute relative error in the computation of 7 ^ with k = 1 for different truncation 
index s; that is, with the approximation 7 ^ — —TlJi. For 4 < s < 6, this error is 
seen to remain essentially constant at 0(10“^). The contribution with /c = 2 is about 
5 orders of magnitude smaller than the k = 1 contribution, so that this additional 
contribution needs to be included for larger index s. The absolute relative error including 
the contribution with k = 2 is shown in the third column; that is, with the approximation 
7 „ ~ — 7y(Ji + J 2 ). The expansion with A: = 3 is about 8 orders of magnitude smaller 
than the k = 1 contribution, so this would only begin to make a significant contribution 
for s > 6 . This problem becomes even more acute for smaller n values, say n = 10, 
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Table 4; Values for 7 „ obtained from (1.1) and (2.13) compared with the exact value. 


n 

Eq. (1.1) 

Eq. (2.13) 

Exact 

In 

10 

+2.105395 X 10-^ 

+2.04713213 X 10“^ 

+2.05332814 

.. X 10-^ 

50 

+1.275493 X 10^ 

+1.26823798 x 10^ 

+1.26823602 

.. X 10^ 

80 

+2.514857 X 10^° 

+2.51633995 x 10^° 

+2.51634410 

.. X 10^° 

100 

-4.259408 X 10^^ 

-4.25340036 x 10^^ 

-4.25340157 

.. X 10^^ 

137 

+3.898740 X lO^'^ 

-7.99377883 x lO^^ 

-7.99522199 

.. X 1027 

200 

-7.060244 X 10®^ 

-6.97465335 x lO^® 

-6.97464971 

.. X 10^5 

500 

-1.165662 X lO^o^ 

-1.16550527 X lO^o^ 

-1.16550527 

. . X 10204 


where higher k values need to be retained. However, the chief interest in the asymptotic 
expansion in ( 2 . 12 ) is for large n, where this problem is of no real concern. 

In Table 4 we show some examples of the asymptotic approximation given in (2.13). 
We compare these with the values produced by the leading approximation (1.1) and the 
exact value of 7 ^ obtained from Mathematica using the command Stielt jesGammafn]. 
It will be observed that for n = 500 the approximation (2.13) yields nine significant 
figures. 

Finally, we remark that the analysis in Section 2 is immediately applicable to the more 
general Stieltjes constants 7 n(cK) appearing in the Laurent expansion for the Hurwitz zeta 
function C('S)Ci) about the point s = 1. These constants are defined by 


where 70 (a) 


C(s,a) 



00 / 'in 

+ M-7n(a) (s - 1) 


—r'(a)/r(a) and 7n(l) = 7n- It is shown in [7, Eq. (2.9)] that 


Cn{a) := 7n(a) 


gU log log a 
a 


00 

k=l 


Then it follows that the expansions in Theorems 1 and 2 are modified only in the argument 
of the trigonometric functions appearing therein, which become na + b — 27ra. Thus, for 
example, from (2.13) we have 

Cnioi) ~ — %=r- {cos (na+6—27ra) f 1 + — + — sin {na+b—2TTa) {— + 1 

( \ n J \n J) 

as n ^ 00 , where the quantities A, B, a, b and the coefficients Cs,ds (s = 1,2) are 
as specified in Theorem 2. The leading approximation agrees with that obtained in [7, 
Eq. (2.4)]. 
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